Prospecção Séries Temporais

G4: Edinaldo de Alencar / Igor Freire / Ramon Araújo / Ricardo Ribeiro

28 de outubro de 2014

Importação da base de dados

Pacote utilizado para leitura em formato de séries temporais:

library(zoo)

Importar as bases de dados utilizando read.zoo:

# A base de dados e os scripts R estão no mesmo diretório (o diretório atual)
setwd(paste("~/Documents/Mestrado/UFPA/Mineração de Dados/data-mining-ppgee/trabalho-2-forecast/", sep=""))
#setwd(paste("./", sep=""));

# Inicialmente, ler a base de dados diários como um data frame (através de read.csv)
dataframe_diario <- read.csv(file = "dataset_diario.csv", sep = ";", dec = ",", header = TRUE)
# Em seguida, converter para uma série temporal (lista indexada pela data)
dataset_diario <- zoo(as.matrix(dataframe_diario[, -1:-2]), as.Date(dataframe_diario[,1], format = "%d/%m/%y"))
# Obs: em "format" usa-se y minusculo, pois a data está no formato dd/mm/yy

# A função read.zoo() abaixo não retorna lista com vetores categóricos e numéricos ao mesmo tempo.
# Isto é, se houver dados categóricos e numéricos na base, os numéricos serão convertidos.
# Por isso, a base com dados diários não foi lida diretamente com read.zoo. A base com dados mensais 
# pode ser lida diretamente com read.zoo()
# Ler o .csv como uma série temporal (indexado pela data)
dataset_mensal <- read.zoo(file = "dataset_mensal.csv", sep = ";", dec = ",", header = TRUE, 
        index = 1, tz = "", FUN = as.yearmon, format = "%m/%Y", drop = FALSE)
# index ->  coluna do arquivo .csv que contém a data
# Obs: em "format" usa-se Y maiúsculo, pois a data está no formato mm/yyyy

Notar que é necessário definir:

  1. O caractere que separa as entradas no arquivo .csv (sep =)
  2. O caractere separador de casas decimais do atributo numérico presente nas bases de dados (dec =)

Informação sobre os dias da semana

Observar que a informação sobre os dias da semana é redundante, pois pode ser obtida através de:

dia <- weekdays(index(dataset_diario))

que tem como resultado, por exemplo, para as 10 primeiras amostras da base:

head(dia, 10)
##  [1] "Terça Feira"   "Quarta Feira"  "Quinta Feira"  "Sexta Feira"  
##  [5] "Sábado"        "Domingo"       "Segunda Feira" "Terça Feira"  
##  [9] "Quarta Feira"  "Quinta Feira"

o qual pode ser confirmado comparando com as 10 primeiras entradas na coluna “dia” da base:

head(dataframe_diario[,2], 10)
##  [1] ter qua qui sex sab dom seg ter qua qui
## Levels: dom qua qui sab seg sex ter

Por isso, na leitura da base em dataset_diario, esta coluna foi ignorada.

Inspeção do Conjunto de Dados

summary(dataset_diario)
##      Index            dataset_diario 
##  Min.   :2002-01-01   Min.   :11130  
##  1st Qu.:2004-01-01   1st Qu.:19468  
##  Median :2005-12-31   Median :22811  
##  Mean   :2005-12-31   Mean   :23158  
##  3rd Qu.:2007-12-31   3rd Qu.:26712  
##  Max.   :2009-12-31   Max.   :34292
summary(dataset_mensal)
##      Index          fluxo       
##  Min.   :1992   Min.   :261544  
##  1st Qu.:1996   1st Qu.:408961  
##  Median :2001   Median :529912  
##  Mean   :2001   Mean   :548440  
##  3rd Qu.:2005   3rd Qu.:662250  
##  Max.   :2010   Max.   :982708

Séries Temporais

Conversão dos dados para o formato de série temporal no R.

Série Mensal:

# Frequency --> número de observações por unidade de tempo 
# define a unidade de tempo (e.g. 12: unidade de tempo = ano)
tsMensal <- ts(dataset_mensal, frequency=12, start=1992)
plot(tsMensal, xlab="Anos", ylab="Fluxo Mensal")

plot of chunk unnamed-chunk-10

Série Diária:

Série temporal diária

tsDiario <- ts(dataset_diario, start=c(2002,1,1), frequency=365)

Série temporal diária multi-sazonal

tsDiarioMSzl <- msts(tsDiario, seasonal.periods=c(7,365.25))
plot(tsDiario, xlab="Anos", ylab="Fluxo")

plot of chunk unnamed-chunk-13

Decomposição das séries temporais

Decompor a série temporal em:

plot(decompose(tsMensal), xlab="Anos")

plot of chunk unnamed-chunk-14

plot(decompose(tsDiarioMSzl), xlab="Anos")

plot of chunk unnamed-chunk-14

  1. Clara tendência (trend) de subida.
  2. Sazonalidade de 1 ano.

Objetivos

Realizar prospecções de curto, médio e longo prazo.

Curto Prazo

Para curto prazo, será utilizada a série temporal com dados diários (tsDiario).

Médio Prazo

Longo Prazo

Para médio e longo prazo, será utilizada a série temporal com dados mensais (tsMensal).

Separação dos dados

Separação das séries temporais em conjuntos de treinamento e de testes.

Curto Prazo

Separação de tsDiario em conjuntos de treinamento e teste.

Conjunto de Treinamento

tsDiarioTrain <- window(tsDiario, end=c(2005,366)) # De Jan 2002 a Dez 2005

Série temporal multi-sazonal do conjunto de treinamento

# Considerar as sazonalidades semanal e anual
tsDiarioTrainMszl <- msts(tsDiarioTrain, seasonal.periods=c(7,365.25))

Conjuntos de Teste

30 dias:
tsDiarioTest30Days <- window(tsDiario, start=c(2006,1), end=c(2006,30)) # De 1/1/2006 a 30/1/2006
45 dias:
# 45 dias a partir De 1/1/2006
tsDiarioTest45Days <- window(tsDiario, start=c(2006,1), end=c(2006,45)) 

Médio/Longo Prazo

4 meses:

Separação de tsMensal em conjuntos de treinamento e teste:

tsMensalTrain <- window(tsMensal, end=c(2005, 12)) # De Jan 1992 a Dez 2005
tsMensalTest4mth <- window(tsMensal, start=2006, end=c(2006,4)) # De Jan 2006 a Mar 2006

6 meses:

tsMensalTest6mth <- window(tsMensal, start=2006, end=c(2006,6)) # De Jan 2006 a Jun 2006

1 ano:

tsMensalTest1yr <- window(tsMensal, start=2006, end=2007) # De Jan 2006 a Jan 2007

2 anos:

tsMensalTest2yr <- window(tsMensal, start=2006, end=(2008)) # De Jan 2006 a Jan 2008

Biblioteca Utilizada

Pacote utilizado (ver (R. J. Hyndman 2014) e (Leek 2014)):

library(forecast)

Estratégia de Avaliação dos Métodos

Apresenta-se a seguir:

Nota:

Saída: objeto forecast, contendo:

Acurácia

Métricas comparativas usuais:

Diferentemente do ME, MAE e RMSE, o medidor MAPE é independente da escala.

MAPE:

\[M = \frac{1}{n} \sum \limits_{t=1}^{n} \left| \frac{A_t - F_t}{A_t} \right|,\]

onde \(A_t\) é o valor verdadeiro, \(F_t\) é o valor predito e \(n\) é o número de amostras temporais preditas.

Como desvantagem, o critério MAPE tem comportamento indefinido quando o valor verdadeiro \(A_t=0\) (divisão por zero).

Métodos simplistas

Funções utilizadas na biblioteca forecast:

Métodos simplistas em prospecções de Curto Prazo

Prospecções através dos métodos simplísticos para um prazo 30 dias.

Série temporal de treinamento

Fluxo diário de Janeiro de 2002 a Dezembro de 2005

# Dados de treinamento
plot(tsDiarioTrain, xlab="Anos", ylab="Fluxo Diário",)
title("Fluxo diário de treinamento: Jan 2002 a Dez 2005")

plot of chunk unnamed-chunk-24

Número de dias para os quais deseja-se prever o fluxo

diasAPrever <- 30

Média

# Média
f_mean <- meanf(tsDiarioTrain, h=diasAPrever)
plot(f_mean, xlab="Anos", ylab="Fluxo Diário", include=365)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-27

# Acurácia
accuracy(f_mean, tsDiarioTest30Days)
##                     ME RMSE  MAE    MPE  MAPE  MASE   ACF1 Theil's U
## Training set 8.412e-13 2461 1997 -1.638 10.47 1.067 0.7655        NA
## Test set     2.813e+03 3246 2923 12.066 12.66 1.562 0.3835     1.848

Naïve

# Naïve
f_naive <- naive(tsDiarioTrain, h=diasAPrever)
plot(f_naive, xlab="Anos", ylab="Fluxo Diário", include=365)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-29

# Acurácia
accuracy(f_naive, tsDiarioTest30Days)
##                    ME RMSE  MAE     MPE  MAPE   MASE    ACF1 Theil's U
## Training set    4.393 1679 1209 -0.3758  6.38 0.6459 -0.1089        NA
## Test set     2058.777 2620 2338  8.6795 10.15 1.2491  0.3835     1.474

Naïve Sazonal

# Naïve Sazonal
f_seasonal_naive <- snaive(tsDiarioTrainMszl, h=diasAPrever)
plot(f_seasonal_naive, xlab="Anos", ylab="Fluxo Diário", include=365)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-31

# Acurácia
accuracy(f_seasonal_naive, tsDiarioTest30Days)
##                ME RMSE  MAE   MPE  MAPE  MASE   ACF1 Theil's U
## Training set 1575 2411 1872 7.265 9.030 1.000 0.1352        NA
## Test set     1464 2280 1757 6.230 7.664 0.939 0.2291     1.403

Drift

# Drift
f_drift <- rwf(tsDiarioTrain, drift = TRUE, h=diasAPrever)
plot(f_drift, xlab="Anos", ylab="Fluxo Diário", include=365)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-33

# Acurácia
accuracy(f_drift, tsDiarioTest30Days)
##                     ME RMSE  MAE     MPE  MAPE   MASE    ACF1 Theil's U
## Training set 1.642e-12 1679 1209 -0.3986 6.382 0.6459 -0.1089        NA
## Test set     1.993e+03 2560 2289  8.3882 9.951 1.2232  0.3755      1.44

Métodos simplistas em prospecções de Médio e Longo Prazo

Prospecções através dos métodos simplísticos para um prazo 4 meses.

Número de meses para os quais deseja-se prever o fluxo

mesesAPrever <- 4

Comparação dos métodos simplísticos

Acuráricas (MAPE)

accuracies
##               Mean Näive N-Sazonal Drift
## Training set 22.82 4.376     7.673 4.408
## Test set     29.64 5.202     6.735 6.139

Métodos simplistas em prospecções de Longo Prazo

Prospecções através dos métodos simplísticos para um prazo 1 ano.

Número de meses para os quais deseja-se prever o fluxo

mesesAPrever <- 12

Comparação dos métodos simplísticos

Acurácias MAPE

accuracies
##               Mean Näive N-Sazonal Drift
## Training set 22.82 4.376     7.673 4.408
## Test set     34.45 5.524     7.163 4.601

Regressão Linear

30 dias

diasAPrever <- 30
reg <- tslm(tsDiarioTrain ~ trend + season)
regfcast <- forecast(reg, h=diasAPrever)
plot(regfcast, xlab="Anos", ylab="Fluxo Diario", include=120)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-42

accuracy(regfcast, tsDiarioTest30Days)
##                      ME RMSE    MAE     MPE  MAPE   MASE   ACF1 Theil's U
## Training set -6.536e-14 1262  984.7 -0.4535 5.253 0.5261 0.2681        NA
## Test set      7.657e+00 1794 1447.1 -0.5431 6.647 0.7732 0.4121    0.9895

1 Ano

mesesAPrever <- 12
reg <- tslm(tsMensalTrain ~ trend + season)
regfcast <- forecast(reg, h=mesesAPrever)
plot(regfcast, xlab="Anos", ylab="Fluxo Diario", include=36)
lines(tsMensalTest1yr, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-43

accuracy(regfcast, tsMensalTest1yr)
##                      ME  RMSE   MAE     MPE  MAPE   MASE   ACF1 Theil's U
## Training set -1.754e-13 26528 21375 -0.3482 4.503 0.5763 0.9182        NA
## Test set      4.926e+04 55298 49262  6.6054 6.605 1.3281 0.6854      1.27

Modelo TBATS

Pesquisar mais.

Método TBATS em prospecções de Curto Prazo

Para 30 dias

diasAPrever <- 30

Modelo, a partir da série multi-sazonal do conjunto de treinamento:

tbats_model <- tbats(tsDiarioTrainMszl)

Predição:

tbats_fc30Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc30Days, include = 120)
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-46

Acurácia:

accuracy(tbats_fc30Days, tsDiarioTest30Days)
##                    ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set    4.494 687.8 393.3 -0.109 2.095 0.2101 0.0006724        NA
## Test set     -210.910 754.8 600.6 -1.086 2.710 0.3209 0.7607279    0.4361

Método TBATS em prospecções de Curto Prazo

Para 45 dias

diasAPrever <- 45

Predição:

tbats_fc45Days <- forecast(tbats_model, h=diasAPrever)
plot(tbats_fc45Days, include = 120)
lines(tsDiarioTest45Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-49

Acurácia:

accuracy(tbats_fc45Days, tsDiarioTest45Days)
##                   ME  RMSE   MAE     MPE  MAPE   MASE      ACF1 Theil's U
## Training set   4.494 687.8 393.3 -0.1090 2.095 0.2101 0.0006724        NA
## Test set     -13.861 722.4 595.0 -0.2168 2.638 0.3179 0.7639143    0.3844

Modelo ARIMA

ARIMA (Autoregressive Integrated Moving Average)

Modelo médias-móveis

Dado no instante \(t\) de uma série temporal é dado pelo valor esperado da variável aleatória acrescida do ruído branco no instante \(t\) e a soma ponderada do ruído branco em instantes passados.

O número de instantes passados considerados é dado pela ordem do modelo.

Médias móveis de ordem 3:

plot(tsMensalTrain, xlab="Anos", ylab="Fluxo Mensal")
movAvg <- ma(tsMensalTrain, order=3)
lines(movAvg, col="red")

plot of chunk unnamed-chunk-51

Auto-regressive - Modelo auto regressivo (AR)

Dado no instante \(t\) de uma série temporal é expresso como a média ponderada dos dados em instantes passados acrescida de ruído branco.

O número de instantes passados considerados é dado pela ordem do modelo.

ARMA

ARMA: autoregressive moving-average - combinação de ambos

ARIMA

ARIMA: ARMA com um passo de diferenciação

Método ARIMA em prospecções de curto prazo

Para 30 dias

diasAPrever <- 30

Modelo:

arima_model_diario <- auto.arima(tsDiarioTrain)

Predição:

fcast_arima_30days <- forecast(arima_model_diario, h=diasAPrever)
plot(fcast_arima_30days, xlab="Anos", ylab="Fluxo Diário", include=120)
lines(tsDiarioTest30Days, col="red")

plot of chunk unnamed-chunk-54

Acurácia:

accuracy(fcast_arima_30days, tsDiarioTest30Days)
##                  ME RMSE    MAE     MPE  MAPE   MASE     ACF1 Theil's U
## Training set  1.700 1198  988.2 -0.3686 5.182 0.5280 0.005709        NA
## Test set     -2.018 1585 1333.9 -0.5408 6.124 0.7127 0.437792    0.8744

Método ARIMA em prospecções de médio prazo

Para 4 meses

mesesAPrever <- 4

Modelo:

arima_model_mensal <- auto.arima(tsMensalTrain)

Predição:

fcast_arima_4mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_4mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest4mth, col="red")

plot of chunk unnamed-chunk-58

Acurácia:

accuracy(fcast_arima_4mth, tsMensalTest4mth)
##                ME  RMSE  MAE     MPE  MAPE   MASE      ACF1 Theil's U
## Training set  145 10236 7610 -0.0156 1.595 0.2052 -0.005073        NA
## Test set     2173 10234 8719  0.2861 1.276 0.2351 -0.511317    0.1839

Para 6 meses

mesesAPrever <- 6

Predição:

fcast_arima_6mth <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_6mth, xlab="Anos", ylab="Fluxo Mensal", include=36)
lines(tsMensalTest6mth, col="red")

plot of chunk unnamed-chunk-61

Acurácia:

accuracy(fcast_arima_6mth, tsMensalTest6mth)
##                ME  RMSE  MAE     MPE  MAPE   MASE      ACF1 Theil's U
## Training set  145 10236 7610 -0.0156 1.595 0.2052 -0.005073        NA
## Test set     3437  9902 8148  0.4719 1.181 0.2197 -0.352187    0.2011

Método ARIMA em prospecções de longo prazo

Para 1 ano

mesesAPrever <- 12

Predição:

fcast_arima_1yr <- forecast(arima_model_mensal, h=mesesAPrever)
plot(fcast_arima_1yr, include=48)
lines(tsMensalTest1yr, col="red")

plot of chunk unnamed-chunk-65

Acurácia:

accuracy(fcast_arima_1yr, tsMensalTest1yr)
##                 ME  RMSE   MAE     MPE  MAPE   MASE      ACF1 Theil's U
## Training set   145 10236  7610 -0.0156 1.595 0.2052 -0.005073        NA
## Test set     22409 31001 24765  2.9237 3.278 0.6677  0.697780    0.7174

Método ARIMA - Conclusões

** Sugiro fazer uma conclusão só, ao final do trabalho, dizendo em quais situações cada método foi mais eficaz. **

Exponential Smoothing

Exponential Smoothing - 2

Ver (R. J. Hyndman and Athanasopoulos 2013). O autor possui um livro exclusivamente sobre o assunto em (R. Hyndman et al. 2008).

ETS

Há, portanto, 30 métodos distintos.

Exponential Smoothing no R

Método Exponential Smoothing em prospecções de Curto Prazo

diasAPrever <- 30

Modelo:

etsDiario <- ets(tsDiarioTrain)
## Warning: I can't handle data with frequency greater than 24. Seasonality
## will be ignored. Try stlf() if you need seasonal forecasts.

Prospecção:

fcastDiario30days <- forecast(etsDiario, h=diasAPrever)

Comparar predição com valores reais do conjunto de teste:

plot(fcastDiario30days, xlab="Anos", ylab="Fluxo Mensal");
lines(tsDiarioTest30Days, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-70

Acurácia:

accuracy(fcastDiario30days, tsDiarioTest30Days)
##                   ME RMSE  MAE     MPE  MAPE   MASE   ACF1 Theil's U
## Training set   80.94 1490 1265 -0.1819 6.702 0.6761 0.3243        NA
## Test set     -169.73 1629 1270 -1.3214 5.917 0.6787 0.3835     0.837

Método Exponential Smoothing em Médio Prazo

Para 6 meses

mesesAPrever <- 6

Modelo:

etsMensal <- ets(tsMensalTrain)

Prospecção:

fcastMensal6mth <- forecast(etsMensal, h=mesesAPrever)

Comparar predição com valores reais do conjunto de teste:

plot(fcastMensal6mth, xlab="Anos", ylab="Fluxo Mensal");
lines(tsMensalTest6mth, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-75

Acurácia:

accuracy(fcastMensal6mth, tsMensalTest6mth)
##                ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set  910  8952  6764 0.1729 1.438 0.1824 -0.005476        NA
## Test set     4511 10353 10213 0.6550 1.495 0.2753  0.152408    0.1901

Método Exponential Smoothing em Longo Prazo

Para 1 ano

mesesAPrever <- 12

Prospecção:

fcastMensal1yr <- forecast(etsMensal, h=mesesAPrever)

Comparar predição com valores reais do conjunto de teste:

plot(fcastMensal1yr);
lines(tsMensalTest1yr, col="red")
legend("topleft",legend=c("Real"),col=2,lty=1)

plot of chunk unnamed-chunk-79

Acurácia:

accuracy(fcastMensal1yr, tsMensalTest1yr)
##                 ME  RMSE   MAE    MPE  MAPE   MASE      ACF1 Theil's U
## Training set   910  8952  6764 0.1729 1.438 0.1824 -0.005476        NA
## Test set     21648 29250 24499 2.8481 3.268 0.6605  0.789771    0.6718

Usando transformação Box-cox

lam <- BoxCox.lambda(tsMensalTrain)
etsMensal_boxcox <- ets(tsMensalTrain, additive=TRUE, lambda=lam)
fcastMensal1yr_boxcox <- forecast(etsMensal_boxcox,  h = mesesAPrever)

plot(fcastMensal1yr_boxcox, include=60)
lines(tsMensalTest1yr, col="red")

plot of chunk unnamed-chunk-81

Acurácia:

accuracy(fcastMensal1yr_boxcox, tsMensalTest1yr)
##                   ME  RMSE   MAE    MPE  MAPE   MASE     ACF1 Theil's U
## Training set   605.1  8996  6818 0.1216 1.456 0.1838 0.008943        NA
## Test set     20121.4 27805 23231 2.6408 3.099 0.6263 0.788149    0.6388

STLF

mesesAPrever <- 12

Dados com sazonalidade removida:

tsMensal.stl <- stl(tsMensalTrain[,1], s.window=12)
# Seasonally adjusted data constructed by removing the seasonal component.
plot(seasadj(tsMensal.stl))

plot of chunk unnamed-chunk-84

stlf_model <- stlf(tsMensalTrain[,1])
stlf_fcast <- forecast(stlf_model, h=mesesAPrever)
plot(stlf_fcast)
lines(tsMensalTest1yr, col="red")

plot of chunk unnamed-chunk-85

Acurácia:

accuracy(stlf_fcast, tsMensalTest1yr)
##                   ME  RMSE   MAE      MPE  MAPE   MASE    ACF1 Theil's U
## Training set  -120.1  7919  6023 -0.02844 1.297 0.1624 0.02359        NA
## Test set     11044.4 22662 18177  1.36140 2.421 0.4901 0.54398    0.5458

Referências

Hyndman, Rob J. 2014. “Forecasting Time Series Using R.” http://robjhyndman.com/talks/MelbourneRUG.pdf.

Hyndman, Rob J, and George Athanasopoulos. 2013. Forecasting: principles and Practice. OTexts.

Hyndman, Rob, Anne B. Koehler, J. Keith Ord, and Ralph D. Snyder. 2008. Forecasting with Exponential Smoothing: The State Space Approach (Springer Series in Statistics). 2008th ed. Springer. http://amazon.com/o/ASIN/3540719164/.

Leek, Jeffrey. 2014. “Practical Machine Learning - Lecture 27 - Forecasting.” Johns Hopkins Bloomberg School of Public Health. https://d396qusza40orc.cloudfront.net/predmachlearn/027forecasting.pdf.